function [f_val_stack, K_store, S_store, theta_mat]=Robust_CS_calc(alpha,a_vec,data,theta_grid_array,Farray,giveMoments,giveWeight,giveVarianceEstim)
%Calculate identification robust confidence sets. theta_grid_array
%is a cell array, where cell i contains the grid of values we want to
%consider to parameter i.  Here I assume we are interested in confidence
%sets for linear combinations of parameters, and the linear combinations of
%interest are given the Farray.  GiveMoments is used to pass the moment and
%derivate functions, while giveWeight is used to pass the desired weighting
%matrix and giveVarianceEstim passes a function used to calculate
%covariance matrices. alpha sets coverage, and a_vec sets the weight to be used
%in linear combination tests (a_vec should be of the same length as the
%third dimension of Farray).

%Set number of simulations to use to critical values
crit_sims=10^5;

%Calculate number of parameters
m=length(theta_grid_array);

%Generate vectors consisting of all parameter combinations consistent with
%the parameter grids
theta_mat=theta_grid_array(1);
theta_mat=cell2mat(theta_mat)';

if m>1
    for n=2:m
        theta_old=theta_mat;
        theta_new=theta_grid_array(n);
        theta_new=cell2mat(theta_new)';
        theta_mat=[kron(ones(length(theta_new),1),theta_old) kron(theta_new,ones(length(theta_old),1))];
    end
end

%Calculate dimension of moments
moments=str2func(giveMoments);
gmat=moments(data,theta_mat(1,:)',0);
k=size(gmat,2);

S_store=zeros(size(theta_mat,1),1);
K_store=zeros(size(theta_mat,1),size(Farray,3));

for n=1:size(theta_mat,1)
    theta=theta_mat(n,:)';
    [S_stat, K_vec]= Robust_stat_calc(data,theta,Farray,giveMoments,giveWeight,giveVarianceEstim);
    S_store(n)=S_stat;
    K_store(n,:)=K_vec';
end

f_val_stack=cell(size(Farray,3),1);
for n=1:size(Farray,3)
    F=Farray(:,:,n);
    F(:,sum(abs(F),1)==0)=[];
    p=size(F,2); %Calculates dimension of tested restriction
    
    a=a_vec(n);
    s = RandStream('mt19937ar','Seed',1);
    RandStream.setGlobalStream(s)
    K_size=sum(randn(crit_sims,p).^2,2);
    J_size=sum(randn(crit_sims,k-p).^2,2);
    crit=quantile((1+a)*K_size+a*J_size,1-alpha);
    init_LC_CS=theta_mat(K_store(:,n)+a*S_store<crit,:);
    f_CS=init_LC_CS*F;
    f_val_stack{n}=f_CS;
end

end

